Pre-visualization of the normalized count data before any differential analysis

Table of content

  1. Sample clustering
    1. With all genes
    2. Without X/Y genes
  2. Check samples and genes with normalized counts
    1. Check for artefact due to FACS sorting
    2. Check for contamination from other cell types that could have escaped the sorting-gating
  3. PCA on the normalized counts
    1. With all factors
    2. SPF vs GF
  4. Gene co-expression analysis
    1. Extract modules of co-expressed genes
    2. Relationship between modules and samples
    3. Genes in modules
      1. Heatmaps with Z-scores
      2. Sinaplots of the Z-scores per groups
      3. Median of counts over ages
    4. Enrichment analysis in modules
    5. Network visualization
    6. Functional gene network

Load

Libraries and functions

Warning message in is.na(x[[i]]):
“is.na() applied to non-(list or vector) of type 'environment'”Warning message in rsqlite_fetch(res@ptr, n = n):
“Don't need to call dbFetch() for statements, only for queries”
==========================================================================
*
*  Package WGCNA 1.63 loaded.
*
*    Important note: It appears that your system supports multi-threading,
*    but it is not enabled within WGCNA in R. 
*    To allow multi-threading within WGCNA with all available cores, use 
*
*          allowWGCNAThreads()
*
*    within R. Use disableWGCNAThreads() to disable threading if necessary.
*    Alternatively, set the following environment variable on your system:
*
*          ALLOW_WGCNA_THREADS=<number_of_processors>
*
*    for example 
*
*          ALLOW_WGCNA_THREADS=4
*
*    To set the environment variable in linux bash shell, type 
*
*           export ALLOW_WGCNA_THREADS=4
*
*     before running R. Other operating systems or shells will
*     have a similar command to achieve the same aim.
*
==========================================================================


Allowing multi-threading with up to 4 threads.
[1] "preparing gene to GO mapping data..."
[1] "preparing IC data..."
[1] "preparing gene to GO mapping data..."
[1] "preparing IC data..."
[1] "preparing gene to GO mapping data..."
[1] "preparing IC data..."

Data

Check samples and specific genes

Check for artefact due to FACS sorting

List of genes to check: Egr1, Jun, Zfp36l1, Malat1, Dusp1, Nr4a1, Fos

Normalized counts

Z-scores

Check for contamination from other cell types that could have escaped the sorting gating

Warning message:
“Setting row names on a tibble is deprecated.”Warning message:
“Setting row names on a tibble is deprecated.”Warning message:
“Setting row names on a tibble is deprecated.”

Remove genes not in list

cell
Scl2a5microglia
Ly6ggranulocyte
Trbc1t
Trbc2t
Igkcb
Ighmb
Mcpt4mast
Cma1mast
NKp46nk
NK1.1nk
NKG2Dnk
Tbetnk

Normalized counts

Z-score

Filter contamination genes and contaminated samples

Genes that could be artefact due to FACS sorting or contamination from other cell types and samples that were potentially contaminated

  1. 19694
  2. 61
  1. 19662
  2. 58
  1. 58
  2. 4

Sample clustering

With all genes

Without X/Y genes

PCA on the normalized counts

With all factors

SPF vs GF

M, 8w

M, 52w

M, 104w

F, 8w

F, 52w

F, 104w

Gene co-expression analysis

Extract modules of co-expressed genes

Filter genes that have a count > 10 in more than 90% of the samples:

  1. 19662
  2. 58
  1. 10254
  2. 58
pickSoftThreshold: will use block size 4363.
 pickSoftThreshold: calculating connectivity for given powers...
   ..working on genes 1 through 4363 of 10254
   ..working on genes 4364 through 8726 of 10254
   ..working on genes 8727 through 10254 of 10254
   Power SFT.R.sq  slope truncated.R.sq  mean.k. median.k. max.k.
1      1    0.119  0.867          0.942 1040.000  1.04e+03 1630.0
2      2    0.349 -1.080          0.939  335.000  3.18e+02  818.0
3      3    0.662 -1.640          0.977  136.000  1.17e+02  494.0
4      4    0.768 -1.920          0.983   64.300  4.87e+01  331.0
5      5    0.805 -2.050          0.979   33.700  2.18e+01  237.0
6      6    0.813 -2.140          0.969   19.100  1.04e+01  177.0
7      7    0.844 -2.130          0.976   11.500  5.21e+00  136.0
8      8    0.842 -2.150          0.976    7.320  2.72e+00  108.0
9      9    0.864 -2.110          0.983    4.840  1.47e+00   87.0
10    10    0.876 -2.070          0.986    3.300  8.05e-01   71.2
11    12    0.897 -1.990          0.991    1.670  2.65e-01   49.4
12    14    0.911 -1.920          0.995    0.920  9.30e-02   35.5
13    16    0.919 -1.870          0.994    0.540  3.45e-02   26.2
14    18    0.930 -1.810          0.996    0.334  1.34e-02   19.7
15    20    0.940 -1.760          0.997    0.216  5.40e-03   15.1
PowerSFT.R.sqslopetruncated.R.sqmean.k.median.k.max.k.
1 0.1187056 0.8673413 0.9420926 1041.51600891.038301e+031628.48756
2 0.3485840 -1.0845046 0.9390350 335.39986303.177673e+02 817.70489
3 0.6619483 -1.6433270 0.9774235 136.26032401.173439e+02 494.24803
4 0.7679699 -1.9178567 0.9833591 64.27162944.866661e+01 331.03946
5 0.8047170 -2.0548036 0.9792207 33.69804042.178180e+01 236.61191
6 0.8126904 -2.1399321 0.9685043 19.12562581.041547e+01 176.81449
7 0.8438263 -2.1326059 0.9762230 11.54569105.213720e+00 136.45578
8 0.8418240 -2.1499432 0.9755992 7.32260802.721281e+00 107.90242
9 0.8635171 -2.1056619 0.9828832 4.83557561.465510e+00 86.95871
10 0.8758574 -2.0665490 0.9856631 3.30240108.047443e-01 71.15384
12 0.8972470 -1.9887865 0.9910725 1.67039882.650155e-01 49.35077
14 0.9109978 -1.9200346 0.9951679 0.91962369.299832e-02 35.46703
16 0.9187904 -1.8687780 0.9939731 0.54035213.448297e-02 26.17220
18 0.9295070 -1.8101820 0.9961623 0.33431691.344888e-02 19.71336
20 0.9398773 -1.7592677 0.9974440 0.21570515.398326e-03 15.09409
Powermean.k.
5 5 33.6980404
6 6 19.1256258
7 7 11.5456910
8 8 7.3226080
9 9 4.8355756
1010 3.3024010
1112 1.6703988
1214 0.9196236
1316 0.5403521
1418 0.3343169
1520 0.2157051

Parameters for network construction

 Calculating module eigengenes block-wise from all genes
   Flagging genes and samples with too many missing values...
    ..step 1
 ..Working on block 1 .
    TOM calculation: adjacency..
    ..will use 4 parallel threads.
     Fraction of slow calculations: 0.000000
    ..connectivity..
    ..matrix multiplication (system BLAS)..
    ..normalization..
    ..done.
   ..saving TOM for block 1 into file norm_genes_TOM-block.1.RData
 ....clustering..
 ....detecting modules..
 ....calculating module eigengenes..
 ....checking kME in modules..
     ..removing 59 genes from module 1 because their KME is too low.
     ..removing 32 genes from module 2 because their KME is too low.
     ..removing 15 genes from module 3 because their KME is too low.
     ..removing 11 genes from module 4 because their KME is too low.
     ..removing 1 genes from module 5 because their KME is too low.
     ..removing 1 genes from module 6 because their KME is too low.
     ..removing 1 genes from module 7 because their KME is too low.
     ..removing 7 genes from module 8 because their KME is too low.
     ..removing 6 genes from module 9 because their KME is too low.
     ..removing 1 genes from module 10 because their KME is too low.
     ..removing 1 genes from module 11 because their KME is too low.
     ..removing 4 genes from module 13 because their KME is too low.
     ..removing 1 genes from module 15 because their KME is too low.
     ..removing 1 genes from module 16 because their KME is too low.
     ..removing 2 genes from module 17 because their KME is too low.
     ..removing 1 genes from module 20 because their KME is too low.
 ..merging modules that are too close..
     mergeCloseModules: Merging modules whose distance is less than 0.35
       Calculating new MEs...

 ME0  ME1  ME2  ME3  ME4  ME5  ME6  ME7  ME8  ME9 ME10 
2353 1555 1244 1161 1080  955  724  453  390  224  115
10
  • 10 modules with between 115 and 1555 genes
  • 2353 genes not assigned to a module
7901
0
'grey'
1
'turquoise'
2
'blue'
3
'brown'
4
'yellow'
5
'green'
6
'red'
7
'black'
8
'pink'
9
'magenta'
10
'purple'

Relationship between modules and samples

pdf: 2

Genes in modules

Percentages of XY genes per modules

[1] "0: 77 - 3.27241818954526%"
[1] "1: 58 - 3.72990353697749%"
[1] "2: 58 - 4.66237942122186%"
[1] "3: 28 - 2.41171403962102%"
[1] "4: 31 - 2.87037037037037%"
[1] "5: 28 - 2.93193717277487%"
[1] "6: 11 - 1.51933701657459%"
[1] "7: 9 - 1.98675496688742%"
[1] "8: 5 - 1.28205128205128%"
[1] "9: 5 - 2.23214285714286%"
[1] "10: 6 - 5.21739130434783%"

Heatmaps with Z-scores

Type / gender / age

Gender / type / age

Age / type / gender

Sinaplots of the Z-scores per groups

The mean of the Z-score over the samples in the group is plot for each gene

Type / Age

Type / Age / Gender

Age / Type / Gender

Median of counts over ages

Enrichment analysis in modules

[1] "0"

Warning message in grep(txdbPattern, installedPackages):
“argument 'pattern' has length > 1 and only the first element will be used”Warning message in pcls(G):
“initial point very close to some inequality constraints”Fetching GO annotations...
For 1495 genes, we could not find any categories. These genes will be excluded.
To force their use, please run with use_genes_without_cat=TRUE (see documentation).
This was the default behavior for version 1.15.1 and earlier.
Calculating the p-values...
'select()' returned 1:1 mapping between keys and columns

        category over_represented_pvalue under_represented_pvalue numDEInCat
9635  GO:0034660            2.866576e-17                        1        114
9531  GO:0034470            1.112522e-14                        1         87
9621  GO:0034641            3.125920e-14                        1        840
11875 GO:0044260            3.501313e-14                        1       1057
11858 GO:0044237            5.217287e-14                        1       1257
3128  GO:0006399            7.536292e-14                        1         57
      numInCat                                         term ontology
9635       405                      ncRNA metabolic process       BP
9531       301                             ncRNA processing       BP
9621      5271 cellular nitrogen compound metabolic process       BP
11875     6828     cellular macromolecule metabolic process       BP
11858     8369                   cellular metabolic process       BP
3128       158                       tRNA metabolic process       BP
[1] "1"

Warning message in grep(txdbPattern, installedPackages):
“argument 'pattern' has length > 1 and only the first element will be used”Warning message in pcls(G):
“initial point very close to some inequality constraints”Fetching GO annotations...
For 1495 genes, we could not find any categories. These genes will be excluded.
To force their use, please run with use_genes_without_cat=TRUE (see documentation).
This was the default behavior for version 1.15.1 and earlier.
Calculating the p-values...
'select()' returned 1:1 mapping between keys and columns

        category over_represented_pvalue under_represented_pvalue numDEInCat
11875 GO:0044260            6.276070e-22                        1        797
17654 GO:0090304            7.151400e-22                        1        546
5709  GO:0016070            1.097553e-18                        1        490
11359 GO:0043170            1.186845e-18                        1        830
2982  GO:0006139            8.844065e-17                        1        564
14130 GO:0051252            9.562379e-16                        1        419
      numInCat                                             term ontology
11875     6828         cellular macromolecule metabolic process       BP
17654     4158                   nucleic acid metabolic process       BP
5709      3717                            RNA metabolic process       BP
11359     7434                  macromolecule metabolic process       BP
2982      4658 nucleobase-containing compound metabolic process       BP
14130     3034              regulation of RNA metabolic process       BP
[1] "2"

Warning message in grep(txdbPattern, installedPackages):
“argument 'pattern' has length > 1 and only the first element will be used”Warning message in pcls(G):
“initial point very close to some inequality constraints”Fetching GO annotations...
For 1495 genes, we could not find any categories. These genes will be excluded.
To force their use, please run with use_genes_without_cat=TRUE (see documentation).
This was the default behavior for version 1.15.1 and earlier.
Calculating the p-values...
'select()' returned 1:1 mapping between keys and columns

        category over_represented_pvalue under_represented_pvalue numDEInCat
5311  GO:0015031            2.123091e-13                        1        190
12253 GO:0045184            2.796775e-13                        1        204
11873 GO:0044257            5.043272e-11                        1         88
11878 GO:0044265            5.384838e-11                        1        104
7362  GO:0030163            1.810119e-10                        1         98
11644 GO:0043632            1.863056e-10                        1         76
      numInCat                                                   term ontology
5311      1603                                      protein transport       BP
12253     1762                  establishment of protein localization       BP
11873      604                     cellular protein catabolic process       BP
11878      769               cellular macromolecule catabolic process       BP
7362       720                              protein catabolic process       BP
11644      497 modification-dependent macromolecule catabolic process       BP
[1] "3"

Warning message in grep(txdbPattern, installedPackages):
“argument 'pattern' has length > 1 and only the first element will be used”Warning message in pcls(G):
“initial point very close to some inequality constraints”Fetching GO annotations...
For 1495 genes, we could not find any categories. These genes will be excluded.
To force their use, please run with use_genes_without_cat=TRUE (see documentation).
This was the default behavior for version 1.15.1 and earlier.
Calculating the p-values...
'select()' returned 1:1 mapping between keys and columns

        category over_represented_pvalue under_represented_pvalue numDEInCat
17654 GO:0090304            7.881888e-24                        1        423
2982  GO:0006139            1.885316e-21                        1        448
12846 GO:0046483            4.245269e-21                        1        453
5709  GO:0016070            6.383775e-21                        1        380
3382  GO:0006725            5.461445e-20                        1        453
9621  GO:0034641            3.882892e-19                        1        477
      numInCat                                             term ontology
17654     4158                   nucleic acid metabolic process       BP
2982      4658 nucleobase-containing compound metabolic process       BP
12846     4751                    heterocycle metabolic process       BP
5709      3717                            RNA metabolic process       BP
3382      4812     cellular aromatic compound metabolic process       BP
9621      5271     cellular nitrogen compound metabolic process       BP
[1] "4"

Warning message in grep(txdbPattern, installedPackages):
“argument 'pattern' has length > 1 and only the first element will be used”Fetching GO annotations...
For 1495 genes, we could not find any categories. These genes will be excluded.
To force their use, please run with use_genes_without_cat=TRUE (see documentation).
This was the default behavior for version 1.15.1 and earlier.
Calculating the p-values...
'select()' returned 1:1 mapping between keys and columns

        category over_represented_pvalue under_represented_pvalue numDEInCat
12052 GO:0044710            7.483744e-25                        1        327
4335  GO:0009056            2.354925e-21                        1        182
12054 GO:0044712            3.393081e-21                        1        104
11858 GO:0044237            2.932534e-20                        1        630
18743 GO:1901575            3.422357e-20                        1        174
3933  GO:0008152            9.440237e-20                        1        678
      numInCat                                term ontology
12052     3276   single-organism metabolic process       BP
4335      1531                   catabolic process       BP
12054      635   single-organism catabolic process       BP
11858     8369          cellular metabolic process       BP
18743     1476 organic substance catabolic process       BP
3933      9196                   metabolic process       BP
[1] "5"

Warning message in grep(txdbPattern, installedPackages):
“argument 'pattern' has length > 1 and only the first element will be used”Fetching GO annotations...
For 1495 genes, we could not find any categories. These genes will be excluded.
To force their use, please run with use_genes_without_cat=TRUE (see documentation).
This was the default behavior for version 1.15.1 and earlier.
Calculating the p-values...
'select()' returned 1:1 mapping between keys and columns

        category over_represented_pvalue under_represented_pvalue numDEInCat
11858 GO:0044237            1.212214e-27                        1        551
3933  GO:0008152            2.325952e-27                        1        596
9621  GO:0034641            3.652663e-24                        1        385
3552  GO:0007005            2.035399e-23                        1        102
3427  GO:0006807            4.243249e-22                        1        395
11859 GO:0044238            7.259044e-21                        1        524
      numInCat                                         term ontology
11858     8369                   cellular metabolic process       BP
3933      9196                            metabolic process       BP
9621      5271 cellular nitrogen compound metabolic process       BP
3552       610                   mitochondrion organization       BP
3427      5576          nitrogen compound metabolic process       BP
11859     8342                    primary metabolic process       BP
[1] "6"

Warning message in grep(txdbPattern, installedPackages):
“argument 'pattern' has length > 1 and only the first element will be used”Warning message in pcls(G):
“initial point very close to some inequality constraints”Fetching GO annotations...
For 1495 genes, we could not find any categories. These genes will be excluded.
To force their use, please run with use_genes_without_cat=TRUE (see documentation).
This was the default behavior for version 1.15.1 and earlier.
Calculating the p-values...
'select()' returned 1:1 mapping between keys and columns

        category over_represented_pvalue under_represented_pvalue numDEInCat
4681  GO:0009987            1.317029e-10                        1        609
13827 GO:0050794            3.533914e-10                        1        456
13823 GO:0050789            6.303041e-09                        1        469
15856 GO:0065007            1.190885e-08                        1        487
5758  GO:0016192            1.204964e-08                        1         89
11858 GO:0044237            1.344211e-08                        1        426
      numInCat                             term ontology
4681     12997                 cellular process       BP
13827     8769   regulation of cellular process       BP
13823     9254 regulation of biological process       BP
15856     9722            biological regulation       BP
5758      1055       vesicle-mediated transport       BP
11858     8369       cellular metabolic process       BP
[1] "7"

Warning message in grep(txdbPattern, installedPackages):
“argument 'pattern' has length > 1 and only the first element will be used”Warning message in pcls(G):
“initial point very close to some inequality constraints”Fetching GO annotations...
For 1495 genes, we could not find any categories. These genes will be excluded.
To force their use, please run with use_genes_without_cat=TRUE (see documentation).
This was the default behavior for version 1.15.1 and earlier.
Calculating the p-values...
'select()' returned 1:1 mapping between keys and columns

        category over_represented_pvalue under_represented_pvalue numDEInCat
3522  GO:0006955            1.223393e-23                        1         82
1045  GO:0002376            1.221234e-20                        1        117
3519  GO:0006952            1.222474e-17                        1         79
12200 GO:0045087            7.751151e-14                        1         44
11388 GO:0043207            3.731855e-13                        1         50
14385 GO:0051707            3.731855e-13                        1         50
      numInCat                                 term ontology
3522       974                      immune response       BP
1045      1958                immune system process       BP
3519      1148                     defense response       BP
12200      499               innate immune response       BP
11388      679 response to external biotic stimulus       BP
14385      679           response to other organism       BP
[1] "8"

Warning message in grep(txdbPattern, installedPackages):
“argument 'pattern' has length > 1 and only the first element will be used”Fetching GO annotations...
For 1495 genes, we could not find any categories. These genes will be excluded.
To force their use, please run with use_genes_without_cat=TRUE (see documentation).
This was the default behavior for version 1.15.1 and earlier.
Calculating the p-values...
'select()' returned 1:1 mapping between keys and columns

        category over_represented_pvalue under_represented_pvalue numDEInCat
9621  GO:0034641            6.572669e-15                        1        179
3427  GO:0006807            1.104907e-14                        1        186
17654 GO:0090304            1.979359e-13                        1        145
2982  GO:0006139            2.619880e-13                        1        158
12846 GO:0046483            3.796785e-13                        1        160
3382  GO:0006725            1.257606e-12                        1        160
      numInCat                                             term ontology
9621      5271     cellular nitrogen compound metabolic process       BP
3427      5576              nitrogen compound metabolic process       BP
17654     4158                   nucleic acid metabolic process       BP
2982      4658 nucleobase-containing compound metabolic process       BP
12846     4751                    heterocycle metabolic process       BP
3382      4812     cellular aromatic compound metabolic process       BP
[1] "9"

Warning message in grep(txdbPattern, installedPackages):
“argument 'pattern' has length > 1 and only the first element will be used”Fetching GO annotations...
For 1495 genes, we could not find any categories. These genes will be excluded.
To force their use, please run with use_genes_without_cat=TRUE (see documentation).
This was the default behavior for version 1.15.1 and earlier.
Calculating the p-values...
'select()' returned 1:1 mapping between keys and columns

[1] category                 over_represented_pvalue  under_represented_pvalue
[4] numDEInCat               numInCat                 term                    
[7] ontology                
<0 rows> (or 0-length row.names)
[1] "10"

Warning message in grep(txdbPattern, installedPackages):
“argument 'pattern' has length > 1 and only the first element will be used”Fetching GO annotations...
For 1495 genes, we could not find any categories. These genes will be excluded.
To force their use, please run with use_genes_without_cat=TRUE (see documentation).
This was the default behavior for version 1.15.1 and earlier.
Calculating the p-values...
'select()' returned 1:1 mapping between keys and columns

        category over_represented_pvalue under_represented_pvalue numDEInCat
3139  GO:0006412            2.008733e-32                        1         44
11288 GO:0043043            6.423539e-32                        1         44
11625 GO:0043604            1.010377e-30                        1         44
3213  GO:0006518            2.855574e-29                        1         44
11624 GO:0043603            1.096834e-28                        1         45
18738 GO:1901566            6.495308e-24                        1         45
      numInCat                                         term ontology
3139       556                                  translation       BP
11288      575                 peptide biosynthetic process       BP
11625      639                   amide biosynthetic process       BP
3213       691                    peptide metabolic process       BP
11624      798             cellular amide metabolic process       BP
18738     1097 organonitrogen compound biosynthetic process       BP

Network visualization

TOM calculation: adjacency..
..will use 4 parallel threads.
 Fraction of slow calculations: 0.000000
..connectivity..
..matrix multiplication (system BLAS)..
..normalization..
..done.
0610007P14Rik0610009B22Rik0610009O20Rik0610010K14RikMettl260610012G03Rik0610030E20Rik0610037L13RikFmc11110002L01Rik⋯Zw10ZwintZxdaZxdbZxdcZyg11bZyxZzef1Zzz3Hikeshi
0610007P14Rik1111001010⋯0100000111
0610009B22Rik1110000010⋯0111001001
0610009O20Rik1110000100⋯1111011011
0610010K14Rik1001111010⋯0010011110
Mettl260001110010⋯1010010000
0610012G03Rik0001110010⋯1010000000